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A biologically motivated model for spatio-temporal coexistence of two competing species is studied 
by mean-field theory and numerical simulations. In d > 2 dimensions the phase diagram displays an 
extended region where both species coexist, bounded by two second-order phase transition lines be- 
longing to the directed percolation universality class. The two transition lines meet in a multicritical 
point, where a non-trivial critical behavior is observed. 
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INTRODUCTION 



Nonequilibrium models with phase transitions into ab- 
sorbing states arise in studies of biological populations, 
chemical reactions, spreading of diseases, and many 
other [1-5]. An absorbing phase is a subspace of states 
that can be reached but not be left by the dynamics. 
In biologically motivated models absorbing phases occur 
whenever a system reaches a configuration or a dynam- 
ical state from where it cannot escape. In many cases 
an absorbing state is characterized by a complete loss of 
activity. It is believed that phase transitions into an ab- 
sorbing state generically belong to the directed percola- 
tion (DP) universality class [6-8]. DP is extremely robust 
and includes even models with several species of parti- 
cles [9-11] and infinitely many absorbing states [12, 13]. 
Exceptions from DP are usually observed in models with 
long-range interactions, quenched randomness, and non- 
conventional symmetries such as a Z2-symmetry or con- 
servation of particle number. 

In models with several species of particles and sev- 
eral absorbing states, there exist in addition to transi- 
tions where one species goes extinct, multicritical points 
where several reduced control parameters vanish simul- 
taneously. Such multicritical points are characterized by 
a hierarchy of order parameter exponents with only one 
of them being the DP exponent, and by a crossover ex- 
ponent (f) [14, 15]. Again, the multicritical behavior is 
different in the presence of symmetries [16]. 

Most DP models have in common that active sites may 
become inactive at a certain rate, i.e. particles can dis- 
appear spontaneously. But even models with a different 
mechanism for particle removal may belong to the DP 
class. For example, in branching and annihilating ran- 
dom walks particles can only annihilate pairwise, leading 
to an algebraic decay in the subcritical phase. Until re- 
cently, it was believed that such models do not show a 
DP-like phase transition in dimensions larger than two. 
However, as has been shown by Canet et al [17], in higher 
dimensions such transitions occur also, but only beyond a 
minimum annihilation rate, so that freshly produced off- 
spring can annihilate with its parent with a sufficiently 



large probability before diffusing away. 

In this article, we consider a two-species model, which 
is inspired by biology. A special feature of this model is 
that one of the two species (denoted as S) cannot vanish 
spontaneously, instead it can only be destroyed by the 
other one (called A). In one dimension it turns out that 
the two species cannot coexist and the phase transition is 
discontinuous for the first species. In higher dimensions 
the two species can coexist in a certain region of the phase 
diagram, giving rise to three different phases. We find 
numerically that all phase boundaries where one species 
becomes extinct, are in the DP class, while the multicrit- 
ical point where the two critical lines meet, shows a more 
complex scaling behavior. Near the multicritical point, 
the critical region becomes very small, and therefore the 
DP behavior can hardly be seen in numerical simulations. 
We apply numerical studies as well as mean-field calcu- 
lations and scaling arguments. In addition to the sta- 
tionary state, we also investigate the dynamical behavior 
starting from a random initial state. Wc find that two 
different scaling forms are needed to describe different 
dynamical regimes. 

The outline of this paper is as follows: In the next 
section, we motivate and define our model. In Section 3 
we formulate a mean field theory and discuss the corre- 
sponding phase diagram. Numerical simulations in one 
and two spatial dimensions are presented in Section 4. 
In Section 5 we suggest phenomenological scaling theory 
for the dynamical and the stationary regimes. Finally, 
we summarize our findings in Section 6. 



DEFINITION AND MOTIVATION OF THE 
MODEL 



The model is defined on a d-dimensional hypercubic 
lattice with L'^ sites. Each lattice site is occupied either 
by species A, or by species S, or it is empty. The system 
evolves random-sequentially according to the following 
dynamic rules: 

• Any site that has a nearest neighbor occupied by 
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TABLE I: Critical percolation thresholds \^ /fJ- of the con- 
tact process of yl-particles and the corresponding critical ex- 
ponents of DP in d+1 dimensions. 

species A becomes itself occupied by species A with 
rate A^i, irrespective of its previous state. 

• An empty site that has a nearest neighbor occupied 
by species S becomes itself occupied by species S 
with rate A5. 

• A site occupied by species A turns into an empty 
site with rate fi. 



This means that the model follows the reaction-diffusion 
scheme 
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TABLE IL Fixed points of the mean-field equations (2)-(3) 
and the corresponding conditions for stability. 



patches of species S can only be destroyed if they are 
invaded by species A. Since the distribution of A is not 
homogeneous but highly correlated, patches occupied by 
S are not destroyed randomly, and it is not obvious that 
the phase transition to extinction of S has to be in the 
universality class of directed percolation. We will see be- 
low that this is nevertheless the case in dimensions larger 
than one, although the DP behavior is difficult to see 
when the density of A is low. 



MEAN FIELD THEORY 
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The model can be considered as a patch occupancy model 
for the coexistence of two consumer species feeding on 
the same resources. The two species cannot coexist in 
the same patch, and therefore each patch is occupied at 
most by one of the two species [18]. Species A has the 
higher fitness and therefore displaces species S when it 
enters a patch occupied by S. However, species A has a 
certain risk of dying out in a patch, for instance because 
it cannot respond well to adverse circumstances, such as 
illness or extreme weather conditions. A and S may stand 
for an asexual and a sexual species, for instance oribatid 
mites in the soil [19], or ostracods in ponds lakes [20], 
where asexual and sexual species are known to have co- 
existed for a long time. Asexual species have a higher 
reproduction rate, but can accumulate deleterious muta- 
tions [21], which reduce their fitness below that of the 
sexual variant. 

In our model the dynamics of species A is not affected 
by the presence of S since empty sites and sites occu- 
pied by S can both be conquered by A. The dynamics of 
species A is therefore a version of DP, and the phase tran- 
sition of species A to extinction is in the DP universal- 
ity class. More specifically, the dynamic rules of the A's 
are equivalent to those of a so-called contact process [22] , 
which is a well-studied lattice model of DP with random- 
sequential updates. Precise estimates of the percolation 
threshold X^ j [i [23-25] and the critical exponents [5] can 
be found in the literature and are summarized in Table I. 

In absence of species A, species 5* spreads and eventu- 
ally occupies the entire system. In presence of species A 



Let us first discuss the mean-field theory of the model. 
The mean-field equations for the densities of species A 
and S are 

PA{t) = \APA{t){^~PA{t)\-iXpA{t), (2) 

Ps(t) - \sPs{i)\y~PAit)~Ps{i)\~\APs{t)pA{t)- 

(3) 

The mean-field approximation can be viewed as the case 
where the A ox S can move to any patch of the lattice 
with a given rate Xa/L'^ or Xs/L"^, and not only to near- 
est neighbors, in the limit of infinite system size L — > 00. 



A. Stationary case 



In order to determine the mean-field phase diagram we 
computed the stationary solutions of Eqs. (2)- (3). Ta- 
ble II lists the fixed points and the corresponding sta- 
tionary densities pA and ps together with the conditions 
for their stability. The fixed point Fi is unstable since As 
and Xa cannot be negative. The other three fixed points 
are stable in different parts of parameter space. 

In the phase diagram the three fixed points correspond 
to three different phases (see Fig. 1). In the subcritical 
phase (marked by S) species A dies out so that species S 
eventually conquers the whole system, approaching the 
stationary density ps = 1- In the coexistence phase 
(marked by A;S) both species can coexist, leading to a 
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FIG. 1: Mean-field phase diagram. In phase S species A 
dies out so that species S conquers the whole system, while 
in phase A species S dies out. In phase A;S both species 
coexist, leading to a non-trivial fluctuating stationary state. 



non-trivial fluctuating stationary state with the densities 
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(4) 
(5) 



Finally, in the phase denoted by A in Fig. 1, the density 
of species A is so high that species S becomes extinct. 
Here Eq. (4) is still vahd while /05 = 0. 

The three phases are separated by two lines of contin- 
uous phase transitions, a vertical one at Xa/ P ~ 1 and a 
curved one in form of a segment of a parabola 



A^ 






1) 



(6) 



Approaching the vertical line from right to left the sta- 
tionary density pA vanishes linearly. The critical expo- 
nent (3 is therefore 1, just as in the standard mean- field 
theory of directed percolation. The same applies to the 
density ps approaching the curved phase transition line 
from the left. This supports our numerical findings (see 
below) that both lines represent DP transitions. 

The two phase transition lines meet at the multicritical 
point Xa/ P = 1; A5 = 0. In order to discuss the behavior 
in the neighborhood of the multicritical point in more 
detail, it is convenient to introduce the parameters 



EA = Xa/p - 1 : 



es = Xs/p, 



(7) 



which measure the reduced distances from the multicriti- 
cal point in horizontal and vertical direction, respectively. 
In terms of these parameters the curved transition line (6) 
is given by 



1(1 + eA) 



(8) 



Close to the multicritical point, where both parameters 
are small, the densities scale as 
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(9) 
(10) 



As we decrease eA from es to 0, we move horizontally 
from the curved to the vertical transition line, and ps 
increases linearly from to 1. Similarly, as we increase 
es from eA to a value much larger than e^, we move 
from the curved transition line vertically upwards, and 
Ps increases from to a value close to 1. Denoting by 



eA and Ac 



es 



(11) 



the horizontal and vertical distances between a given 
point (eA.es) and the curved transition line, we have 
for A^ <eA,A5<es 



Ps 



Aa 
es 



A5 
eA ' 



(12) 



Therefore, the linear decrease of ps does not depend on 
the direction in which the curved phase transition line 
is approached. However, the slope increases and finally 
diverges as we approach the multicritical point. This 
indicates that the transition with respect to S at the 
multicritical point may no longer belong to the DP class. 



B. Time-dependent mean field solution 

In order to understand the dynamical properties, let us 
now solve the time-dependent mean field equations (2)- 
(3) explicitely. Since species A evolves independently ac- 
cording to the rules of DP Eq. (2) is autonomous. With 
fA = Xa/ p — 1 and es ~ Xs/p the solution (up to a shift 
in time) reads 



PA{t) 



eA 



(1 + eA) [1 - c^p{-eApt)] 



(13) 



Close to criticality in the active phase < ea "C 1 this 
density first decays as a power law 
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(14) 



until it crosses over to the stationary value given in 
Eq. (9) at some typical crossover time tc ~ -^—■ 

Inserting the solution (13) into the second mean field 
equation (3), it is in principle possible to compute the 
time-dependent density ps{t)- 

Let us first consider the case e^ = (the vertical line 
in Fig. 1), where the DP process of species A is critical. 




FIG. 2: Temporal evolution of the density psit) according to 
the mean-field solution (16) for eA = 0, es ~ 0.1, fJ. ~ 1 and 
c = 10^. The three regimes are explained in the text. 



Inserting Eq. (14) into Eq. 
equation reduces to 



(3); the second mean field 
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At the niulticritical point eg = e^i = this differential 
equation further reduces to ps{t) = ~t~^Ps{t) so that 
the density ps{t) ~ 1/i decays in the same way as pAit). 
For €s > 0, however, we find the formal solution 



Psit) = 



exp{espt) 



(esfit) 



l+es 
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esfJ-t 



e^ z 



dz 



(16) 



where C is an integration constant. To get an impression 
about the behavior for es > 0, we plotted this function in 
Fig. 2. As can be seen, there are three different dynamical 
regimes enumerated by I, II, and III. In the first regime 
the density ps{t) decreases algebraically as l/<^+'^', fol- 
lowed by a quick increase in the second regime until the 
density saturates at 1 in regime III. 

The observed behavior can be interpreted as follows. 
Starting with a homogeneously distributed mixture of 
species A and S, the displacement of 5 by ^ initially 
dominates the dynamics. In regime II the spreading of S 
to empty patches, as expressed by the linear term Xspsit) 
on the r.h.s. of Eq. (3), becomes relevant so that the S"- 
population begins to grow exponentially in the voids of 
the critical A-process. The growth of the S'-population 
continues until all accessible voids are filled and the den- 
sity ps{t) reaches its stationary maximum in regime III. 

Next, let us consider the dynamics on the curved phase 
transition line close to the multicritical point. Using 
Eqs. (8) and (13), expanding es and pA in powers of 
eA and keeping only terms up to the third order in the 
densities and/or the control parameter, Eq. (3) turns into 
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FIG. 3: Space-time plot of the 1-1- 1-dimensional model. White 
and black sites correspond to species S and A respectively, 
and grey sites correspond to empty patches. 



The solution is 
psit) = 



[1 - cxjp{-eApt)] 



-l-eA 



C + /i '' (l-e-^)- 



'^ dz 



(18) 



with an integration constant C. This is a monotonous de- 
cay of the density ps- For short times, it is approximated 
by PS '^ 1/^7 as we have obtained before. For large times, 
the density ps approaches the behavior ps — 1/eAfJ-t, i.e., 
the decay is delayed. This large-time behavior is inde- 
pendent of the initial value of pA, and it is therefore also 
obtained when would have started with the stationary 
solution for pA- 



NUMERICAL SIMULATIONS 



Simulations in 1 + 1 dimensions 



The 1+1-dimcnsional case is special in so far as species 
A and species S cannot coexist since the world lines of 
A and S cannot cross each other. Hence a species that 
persists forever leaves no room for a world line of the 
other species. Fig. 3 shows a space-time plot oi A, S and 
empty sites in the parameter range where A persists. As 
can be seen, starting from a random initial state, species 
S eventually dies out. 

The stationary densities for A and S as functions of 
the control parameter Xa/ P were measured in numerical 
simluations using a system size of L = 50 000 sites. As 
shown in Fig. 4, species S undergoes a first-order phase 
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FIG. 4: Stationary density as a function of the control pa- 
rameter X/n in one dimension (here A = Aa = As). 
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FIG. 5: Numerically determined phase diagram in 2+1 di- 
mensions. The dotted hues denote the size of the critical 
regions for variation of Ayi//i(dotted line) and As/^ (dashed 
line). 



transition from unit density to zero density at the critical 
point A^^ of species A. 



B. Stationary properties in 2+1 dimensions 



In higher spatial dimensions both species can coex- 
ist and one obtains a phase diagram with two continu- 
ous transition lines, which is qualitatively similar to that 
of the mean-field approximation. The numerically de- 
termined phase diagram in 2+1 dimensions is shown in 
Fig. 5. The critical thresholds obtained by varying Xa/ IJ- 
agrees with those obtained by varying Xs/ ^J■, which con- 
firms the correctness of the results. The multicritical 
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FIG. 6: Logarithmic plot of the stationary density ps for 
Ag/^ — 60 as Xa/i-i. is varied. A power-law fit in the interval 
IAa/m - (Aa/a*)!* <= 0.05 yields the estimate /3s = 0.583(3). 
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where the evolution of species A is critical and branching 
of S is not allowed. In fact, even for very small A5 > 
species S should be able to survive as species A becomes 
extinct, eventually approaching ps = I. Therefore, as 
in mean-field theory, the multicritical point is located 
exactly at A^^^ = 0. 

In order to determine the critical exponents in the sta- 
tionary state, we performed simulations at several points 
in the phase diagram close to the phase boundaries. The 
lattice size was L^ = 300 x 300 sites. After 200 000 time 
steps the densities of A and S reached a stationary value. 
The fluctuations were within 5% of this value. To deter- 
mine the mean densities, the density was averaged over 
an additional 100 000 time steps. 

For the extinction of species A at the vertical tran- 
sition line, the expected directed percolation behavior 

PA ~ {Xa - X^2')'^^ with Pa = Pdp ^ 0.583 is clearly 
seen in the simulations. However, the evaluation of the 
critical exponent Ps in the vicinity of the curved phase 
transition line turned out to be more difficult. For ex- 
ample, close to the multicritical point the estimates were 
found to depend on the direction in which the line is ap- 
proached, which is known to be impossible in standard 
critical phenomena. In fact, reliable estimates could only 
be obtained far away from the multicritical point. Choos- 
ing Xs/fJ. = 60, where the density of species A is high, 
and estimating the exponent only in a limited range of 
the data close to criticality (see Fig. 6), we find numerical 
evidence that the transition along the curved line belongs 
to the DP universality class. 

Approaching the multicritical point deviations from 
the expected power law behavior become more pro- 
nounced. This is demonstrated in Fig. 7, where we mea- 
sured Ps for Xa/ P = 1, approaching the curved transi- 
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FIG. 7: Logarithmic plot of ps, including its fluctuations for 
Aa/m = 1- The power-law fit gives the exponent 13 = 0.584(9). 



tion line vertically from above. To estimate the size of 
the critical region, the data points for ps were plotted 
including error bars. Then, a power-law was fitted to the 
first k data points, with k being small enough to be close 
to the scaling regime. The measured critical exponent 
was compatible with (3s = 0.585 ±0.05 everywhere. This 
confirms that the entire phase transition line (the solid 
curved line in Fig. 5) except for the multicritical point 
belongs to the DP universality class. 

In order to quantify how the critical behavior is ap- 
proached we defined critical regions whose size is deter- 
mined by the requirement that the error bar of the last 
data point in Fig. 7 still touches the power-law fit. The 
size of these regions with respect to variations of A^ and 
\s are indicated in Fig. 5 by dashed lines. Although their 
size is arbitrary in the sense that it depends on the cho- 
sen accuracy of the computer simulations, it illustrates 
why estimates far away from the multicritical point are 
more reliable. Moreoever, it is interesting to note that 
the critical region for variation of A^//x is much larger 
than the one for Xs/p-- 

As we approach the multicritical point the critical re- 
gions become narrower for both variations. In the imme- 
diate vicinity of the multicritical point (Xa/ fJ- < 0.5 and 
Xs/p. < 1), a reliable estimation of the critical exponent 
(3s becomes extremely difficult due to huge fiuctuations 
of Ps- Here we have an almost instant increase in the 
density of S as we vary the control parameters. Below, 
we will determine the shape of the curved phase tran- 
sition line close to the multicritical point with a better 
precision by using dynamical simulations. 



C. Dynamical properties in 2+1 dimensions 



In order to compare the time dependence of ps{t) in 
2+1 dimensions with the mean field prediction shown in 




FIG. 8: Decay of the densities pA{t) and ps{t) at the mul- 
ticritical point in 2+1 dimensions. The dashed line indicate 
the slopes —5 and —9 (see text). 



Fig. 2, we first performed numerical simulations along the 
vertical line Xa = X]^ , where the DP process of species A 
is critical. Starting with a homogeneous random mixture 
of both species, we monitor the evolution of the densities 
PAit) and ps{t) for various values of Xs/p- 

The results of simulations at the multicritical point 
Xa — X^ ,Xs = arc shown in Fig. 8. As expected, the 



density of species A decays as 
PA{t)^t- 



(20) 



where S ~ P/i^W ~ 0.451 is the usual decay exponent of 
DP in two dimensions. Similarly, the density of species 5" 
is found to decay algebraically as ps{t) ^ t~^ with an 
exponent » 1.6. 

The value of the exponent 9 can be explained as fol- 
lows: Initially each site is randomly occupied by A or 5 
with equal probability. As time evolves the critical DP- 
process of species A gradually removes species 5, leading 
to a monotonous decrease of ps{t)- As species S does not 
create offspring at the multicritical point, all sites occu- 
pied by S at time t are precisely those sites that have 
never been occupied before by species A. Hence ps(t) 
is the probability that a given site has never been vis- 
ited by the DP process of species A. In the literature 
this probability is known as the local persistence prob- 
ability [26] Pi{t) of directed percolation. This quantity 
was shown to decay algebraically as Pi{t) ~ t~^ , where 
6 is the so-called local persistence exponent which seems 
to be independent of the other DP exponents. In one 
spatial dimension its value was estimated by 9 ~ 1.50(1), 
while in two dimensions we find a slightly higher value 
6 = 1.58(3). This observation allows us to conclude that 
Psit) decays at the multicritical point as 



psit) ^ t- 



{eA = es = 0) 



(21) 



Increasing As (while keeping the ^-process critical) we 
allow species S to expand in the voids of the ^-process. 
In this case the algebraic decay (21) is expected to cross 



over to a growth of surviving S'-domains. This density 
continues to increase until the entire empty space is con- 
quered so that PS approaches the value 1. Our simulation 
results for this process are shown in the top left panel of 
Figure 10. Obviously, the temporal behavior oi ps{t) re- 
sembles the mean field solution in Fig. 2, where three 
different dynamical regimes were identified. However, as 
an important difference we note that in regime If the den- 
sity ps increases algebraically as f^ in finite-dimensional 
systems whereas in the mean-field limit the increase is 
exponential. Consequently, the two crossover time scales 
between the regimes are expected to scale differently. 

Next, we performed numerical simulations on the 
curved phase transition line. The results arc shown in the 
top right panel of Figure 10. For short times, ps decays 
again as t~^ , i.e., the difference between the two critical 
lines is not yet felt. For large times, when pA has become 
stationary, ps decays as t~^ , i.e., it behaves as a directed 
percolation process at its critical point. Between these 
two regimes, there is an increase in the density ps, which 
is not seen in the mean-field theory. As we will discuss 
further below, the reason for this is that the characteristic 
time scales for the cutoff of the initial decay and the onset 
of the DP decay scale differently with 65, while there is no 
distinction between the two scales in mean-field theory. 
Therefore, in two- or three-dimensional systems the den- 
sity Ps increases in the intermediate regime, where the 
spreading to empty patches has become relevant, while 
the critical DP behavior, which requires stationary oi pA, 
is not yet established. In analogy to the behavior on the 
vertical transition line, we denote the three regimes of 
initial decay, intermediate increase and ultimate decay 
with regime I, II, and III. 



D. Functional form of the curved transition line 



In the mean-field case the curved transition line was 
shown to be a segment of a parabola which terminates at 
the multicritical point with slope 1 (see Eq. (6)). In low 
dimensions it is therefore natural to assume a power-law 
of the form 



where again 



4«(eAf : 



^W 



CA = (Aa - Ay)//i, es = As/yu 



(22) 



(23) 



are the reduced parameters and y is an unknown ex- 
ponent with the mean- field value y^^ — 1. As in the 
mean-field case, this power law may be superposed by 
next-leading algebraic corrections. 

In order to determine the value of the unknown expo- 
nent y in 2-1-1 dimensions numerically, we adjusted the 
parameters e^ and eg in such a way that ps{t) decays 
eventually in the same way as in a DP process. However, 




FIG. 9: Left panel: Critical line for the extinction of species S 
for small e\ and eg. Right panel; Corresponding local slopes 
for a visual extrapolation. The inset suggests an asymptotic 
value y ~ 0.7. 



the determination of y turns out to be extremely difficult. 
In particular, the estimates seem to drift to smaller val- 
ues as we approach the multicritical point. Extrapolating 
the local slopes (see Fig. 9) we arrive at the conclusion 
that 2/ < 1, probably y w 0.7(2). 



SCALING PROPERTIES 



In this section, we interpret the observed results in 
terms of phenomenological scaling arguments. It turns 
out that the presence of different scales makes it impos- 
sible to formulate a unified scaling theory for the whole 
dynamical range, instead we find two different scaling 
forms for regimes I, II and II,III, respectively with a com- 
mon overlap in region II. 

Starting point is the DP-process of species A which 
decouples from the dynamics of species S. Close to crit- 
icality this DP process is known to be invariant under 
scale transformations 



EA ^ Ae^ , PA ^ A/^ 



PA- 



t^A-'H, 



(24) 



where A > is a scaling factor while f3 and i^n arc DP 
exponents listed in Table I. As usual scale invariance 
implies that the density is a homogeneous function, i.e.. 



PA{eA,t) ^ A-^PA{AeA,A-''H) ^ 



(25) 



Setting A = t^''^\\ and suppressing non-universal metric 
factors one obtains the usual scaling form 



PAit)=t-'FA{eAt'^''^^), 



(26) 



where Fa is a universal scaling function. In order to 
derive a similar scaling form for the dynamics of species S 
in the vicinity of the multicritical point, it is natural to 
postulate analogous scaling properties for £5 and ps, i.e. 



PS 



A^ 



'PS, 



es 



A^^es 



(27) 



with certain exponents xi and X2. Scale invariance then 
implies that 

PsieA, es, t) = A--1 ps{AeA, A^^es, A-'H). (28) 



Choosing again A = t^/'^ii we arc led to the scahng form 
ps{eA,es, t) - i--i/''ii Fs{eA t'/"^^ , es t^'^"^^ ) (29) 

where Fs is another scahng function. However, as we 
will see below, this simple scaling form is not capable to 
describe the whole dynamical evolution of ps , instead it 
can be used only partially either in the regimes I,II or 
in II, III, in each case with a different set of exponents 
Xi and X2- Furthermore, sufficiently close to the curved 
critical line, the distance from the critical line becomes a 
relevant parameter, and we will need a scaling form that 
depends on two variables. 



A. Crossover between the regimes I and II 



In regimes I and II, the decay of pA is yet far from 
becoming stationary, and the decay of ps should be de- 
scribed by the same expression (29) everywhere near the 
multicritical point, and in particular on the two critical 
lines. 

In order to determine the exponents xi and X2 in 
regimes I and II let us first consider the multicritical 
point EA = £s = 0. Here the scaling form (29) re- 
duces to ps{t) ~ t^^i/'^ii Fs(0, 0) which - compared with 
Eq. (21) - gives 



Xl 



Ov« 



(30) 



Moving up the vertical phase transition line, the DP pro- 
cess of species A remains critical but species S recovers 
after some time. It is natural to expect that the crossover 
from decay (regime I) to subsequent growth (regime II) 
takes place at a typical time scale 



f. 



(1) 



1/es 



(31) 



which is inversely proportional to the rate of spreading 
of species S. Scale invariance of the second argument of 
Eq. (29) then implies that 



X2 = v\\. 



Therefore, we arrive at the scaling form 



-0 T?Wl 



psieA.es^t) ^t-" Fl^>{eAt'^''\ est) 



(32) 



(33) 



This scaling form can be verified along the vertical line 
by plotting ps{t)eg versus tes for various values of eg and 
checking for a data collapse. As shown in the second row 
of Fig. 10, this data collapse works nicely at the crossover 
between the dynamical regime I to II (initial decay and 
subsequent increase) while it clearly fails in regime III. 

The first argument of Fg is in fact irrelevant, since 
the density ps does not depend on it in regimes I and II. 
However, it will be relevant in regime HI, where pA be- 
comes stationary and determines the long-term behavior 
of PS- 
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FIG. 10: Simulations in 2+1 dimensions. The left column 
shows results obtained along the vertical transition line ea = 
for as = 0.0001, 0.0002, . . . , 0.1024 while the right column dis- 
plays similar results taken at the curved phase transition line 
for e\ = 0.0016, 0.0032, . . . , 0.2048. Data sets far away from 
the multicritical point (es > 0.0512), where scaling is poor, 
are represented as dotted lines. 

The panels are arranged in three rows. The first row shows 
the densities ps(t) and PA{t), the latter vertically shifted. In 
the second row the scaling form (33) is used to collapse the 
curves in the dynamical regimes I and II (see text). Simi- 
larly, the third row shows data collapses based on the scaling 
form (38) which is valid in the regimes II and III. The respec- 
tive crossovers are indicated by vertical dashed lines. 



B. Crossover between the regimes II and III 



The failure of a data collapse in regime III according 
to Eq. (33) suggests the presence of two different time 
scales in the system, namely, a crossover time scale 



(1) 



(34) 



from where on offspring production becomes relevant, 
and a second time scale tc from where on species S feels 
the stationary behavior of species A which either induces 
a critical DP-like decay of species S towards zero on the 
curved transition line or to an asymptotically stationary 
value of Ps in the coexistence region and on the vertical 
line. To determine this second crossover time let us first 



consider the vertical phase transition Une, where species 
A evolves as a critical contact process. After the initial 
persistence-like decay the remaining S particles are the 
seeds of d-dimensional spheres growing linearly by off- 
spring production essentially unchallenged by species A. 
Therefore, in finite dimensions d the density ps increases 
as f^. On the vertical phase transition line this growth 
continues until the density of species S saturates at the 
value p5 w 1, which defines the second crossover time. 
Thus, matching initial decay and subsequent growth this 
second crossover time has to scale as 






(35) 



Therefore, in finite dimensions the two crossover times 

tc and tc scale differently with respect to eg which 
explains the necessity of two separate scaling forms. 

The exponents xi and X2 for the crossover from regime 
II to regime III can be determined as follows. On the one 
hand, the density ps saturates at the value I for t ^ oo 
on the vertical phase transition line, independent of eg, 
implying that 



x^ =0. 



(36) 



On the other hand, scale invariance of the argument of 
the scaling function in Eq. (29) requires that 



X2 



11 



i + e/d' 

Therefore, the resulting scaling form reads: 



('2'! 



A/^w 



est 



i/{i+e/d) 



(37) 



(38) 



This scaling form can be verified numerically by a data 
collapse along the vertical transition line, as demon- 
strated in Fig. lOe. In Fig. lOf we demonstrate that 
this data collapse works not only at the vertical transi- 
tion line but also along the curved phase transition line. 
Since the DP-like decay of ps is expected to set in when 
PA becomes stationary, i.e., at time i ~ e^ , it follows 
that the two time scales, namely, the DP correlation time 

^11 ~ e^ " and the second crossover time tc ~ £5 
become identical on the curved critical line. More specif- 
ically, the phase transition line is expected to be charac- 
terized by a constant value of the scale-invariant combi- 



nation e^ /e 



-(l+9/d) 



hence 



e*S^{e*AY: 



y = X2 



i + e/d 



(39) 



Inserting the numerical estimates v\\ = 1.295(6) and 6 = 
1.58(3) we obtain 



y- 0.72(1) 



(40) 



in fair agreement with the numerically extrapolated value 
in Sect. 4D. 
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FIG. 11: Left panel: Same data as in Fig. lOf, but with time 
measured in units of ^y ~ e^ on the horizontal axis. Right 
panel: Simulation with the same parameters, but with a sta- 
tionary initial distribution of A particles. 



For large times, the expression (38) for ps must ap- 
proach an asymptotic value that is independent of time. 
We therefore obtain 



limFf (zi,Z2) 



Z2 






(41) 



This implies that lines of constant density in the es-^A- 
plane are given by constant ratios es / Cj^- 

As we have mentioned, the scaling relation (39) implies 
that the DP-like decay of species S sets in only after its 
density has increased to a considerable value. At this 
time species S has lost the memory of its past, and its 
dynamics is essentially determined by species A. There- 
fore, the time evolution ps{t) is from then on identical to 
the one obtained starting the simulation with a station- 
ary Pa and with a finite and random occupation with 5*. 
This conclusion is supported by comparing the two simu- 
lations, as shown in Fig. 11. The range of the axes is the 
same in both panels, and one can clearly see that for suf- 
ficiently large times the two types of curves agree. This 
agreement of the large-time behavior for the two different 
types of initial conditions occurs also in the mean-field 
theory, as mentioned at the end of Section 3. 



Scaling near the curved transition line in 
regime III 



Finally let us study the asymptotic critical behavior 
of ps in the vicinity of the curved phase transition line. 
Close to the multicritical point we have shown that this 
line can be parametrized by 



a(4) 



(42) 



where a is non-universal metric factor. On the curved 
line the density ps{t) decays asymptotically as i~ , i.e., 
as in a critical DP process. Slightly above in the co- 
existence phase, ps{t) will eventually saturate at a value 
proportional to (A5)'', where A5 = es — f*s denotes again 
the vertical distance from the line. This means that in 
addition to the previously discussed invariance under the 
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transformation 



EA 



KeA, es~^A«es, t ^ A-^H 



(43) 



the system simultaneously has to be invariant under the 
DP-hke scale transformation 



Ac 



f7A. 



t^n-'^H., PS 



n^ 



PS 



(44) 



with another scale factor fl independent of A. The second 
type of scale invariance further constraints the form of 
Ps(eA,es,i), effectively reducing the number of indepen- 
dent arguments in the corresponding scaling form to 1. 
To see this it is convenient to switch from the parameters 
(e^,es) to the parameters (ca^Ss), where 



<5s(eA,es) 



K^a) 



eg 



1 



(45) 



is the normalized distance from the critical line. In terms 
of these new parameters, the scaling form (38) is ex- 
pressed equivalently as 



ps{eA,Ss,t) = GsieAt'^-'^^,Ss). 



(46) 



where Gs{zi,Z2) = Fg (zi, az\{l + Z2))- Obviously, the 
parameter Ss is invariant under the first type of scale 
transformation (43) while it scales as Sg — > QSs under 
the second type of scale transformations (44) . Applying 
the latter to the scaling form (46) yields 

Ps{eA,Ss, t) = n-/' Gs{n-hA i'/"" , n6s) . (47) 

Setting ^ ~ 5^ one obtains 

Ps{tA, Ss,t) = Sl Gs{Ss eA t^'"'' )■ (48) 

In fact, the argument 5s £a t^'^''^ is a scale-invariant ratio 
under both types of scale transformations, proving that 
the off-critical DP-process of species A and the almost- 
critical DP-process of species S can be described in terms 
of a unified scaling theory. We emphasize that this scal- 
ing form is valid only in the temporal regime III close to 
the curved phase transition line, where 5s *C 1. 

The scaling form (48) reproduces the expected DP be- 
havior. Initially the system does not yet feel the influence 
of 5s, hence Gs{z) ~ z~'^ for small z. Later, for 5s > 0, 
the initial DP-like decay of ps crosses over to a stationary 
value, meaning that for large arguments Gs{z) becomes 
constant. Hence ps evolves as 



Ps{iA,5s,t) 



(t/Cii)-* for 1 « i/^ii « (5/" , 



5!^ 
"s 



for t/e||»'^, 



s ' 



(49) 

where ^|| ~ e^ is the stationary correlation length of 
species A. 

Exactly on the curved transition line, we have 5s = 0, 
and Ps therefore decays forever. The unit time is set 



by 5|| , which is the life time of the largest voids of the 
A-process. Close to the curved transition line, a decay 
of Ps with t~ occurs only if the time scale for satura- 
tion, ^\\5g is larger than the crossover time tc , imply- 
ing A5 ^ es, or, equivalently, i5s ^ 1. Otherwise, as 
we have already mentioned, the scaling form (48) is not 
valid, and saturation sets in directly after the growth in 
regime II, without DP-like decay of ps- 

Finally, let us compare the results once more to mean- 
field theory. In mean-field theory, the two crossover times 
scale in the same way, making the distinction between 
two different scaling forms unnecessary and removing 
the intermediate growth in regime II close to the crit- 
ical line, where the DP-like decay of ps becomes visi- 
ble. The values of the exponents are in mean-field the- 
ory P = y = 9 = vu =1. The behavior in regime III 
near the critical line that we found in mean-field theory 
resembles qualitatively the one described in this section. 
In particular, the proportionality factor between ps £md 
As diverges at the multicritical point, and so does the 
factor in front of the i~* decay. 



CONCLUSIONS 



In this paper we have introduced a simple model of 
two competing species A and S feeding on the same re- 
sources and living in "patches" represented by lattice 
sites. Species A is characterized by fast reproduction 
and can therefore displace species S and evolves inde- 
pendently according to the dynamic rules of a contact 
process. Hence it belongs to the universality class of di- 
rected percolation. Special S has a lower reproduction 
rate and is therefore restricted to live in those patches 
where A is absent. In contrast to species A, which can 
die out in a patch, patches occupied by species S do not 
become empty spontaneously. 

In one spatial dimension, the phase diagram comprises 
two different phases where one of the species goes ex- 
tinct, separated by a first-order transition line. In higher 
dimensions an additional mixed phase emerges, bounded 
by two continuous transition lines. This mixed phase 
is characterized by a non-trivial coexistence of the two 
species in the stationary state. [27] 

The curved phase transition line (except for the mul- 
ticritical point) is found to belong to the DP univer- 
sality class. This result may be surprising because S- 
individuals do not vanish spontaneously. On the other 
hand the result is also plausible since species S 'perco- 
lates' in the voids of the supercritical DP of species A. 
Since the spatio-temporal arrangement of species A in- 
volves only short-range correlations these voids are un- 
correlated on large scales, which effectively leads to a DP 
transition. 

At the multicritical point the correlation lengths of 
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the A-process diverges, leading to a non-trivial scaling 
behavior. In this case, starting with a mixed random 
initial conditions, 5-occupied patches mark those lattice 
sites that have been never visited by a species A be- 
fore. Therefore, the density ps decays in the same way 
as the so-called local persistence probability in DP. For 
As > 0, however, this decay crosses over to an increase 
of PS when the reproduction of S becomes visible. This 
increase continues until a stationary state is reached in 



which the two species coexist if the parameters are not 
too close to the curved phase transition line. Close to 
the curved phase transition line, ps decreases again after 
times long enough that pA has become stationary and 
saturates eventually at a small nonzero value depend- 
ing on the distance to the critical line. We succeeded in 
introducing scaling variables and scaling functions that 
characterize the behavior of ps in the different regimes. 
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